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Probing cosmic star formation up to z = 9.4 with GRBs 
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ABSTRACT 

We propose a novel approach, based on Principal Components Analysis, to the use 
of Gamma-Ray Bursts (GRBs) as probes of cosmic star formation history (SFH) up 
to very high rcdshifts. The main advantage of such approach is to avoid the necessity 
of assuming an ad hoc parameterization of the SFH. We first validate the method 
by reconstructing a known SFH from Monte Carlo-generated mock data. We then 
apply the method to the most recent Swift data of GRBs with known redshift and 
compare it against the SFH obtained by independent methods. The main conclusion 
is that the level of star formation activity at z « 9.4 could have been already as high 
as the present-day one (« O.OlA/ yr™ 1 Mpc -3 ). This is a factor 3-5 times higher 
than deduced from high-z galaxy searches through drop-out techniques. If true, this 
might alleviate the long-standing problem of a photon-starving reionization; it might 
also indicate that galaxies accounting for most of the star formation activity at high 
redshift go undetected by even the most deep searches. 
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1 INTRODUCTION 

The cosmic star formation history (SFH) is an important 
test for galaxy formation mod els. Experimentally o ur kn owl- 
edge of the SFH comes from iHopkins fc Beacoml (|2006l ) up 
to z m 6 and from observations of color-selected Lyma n 
Break Galaxies ( Mannucci ct al 1l2007l ;l Bouwens et al 
Lya Emitters |Ota et all 20081 ) . UV+IR measurements 
jReddv et alj 2008 ) . and GRB observations (|Charv et al.l 
120071 ; lYiikselet all [2001 ; IWang fc Pail 120091 ) at higher z 



(hereafter, these will be refereed to as H2006, M2007, B2008, 
O2008, R2008, C2007, Y2008 and W2009, respectively). 
However, direct high-z measurements constitute an extreme 
challenge even for the most powerful telescopes and remain 
sparse. 

Due to their high luminosity GRBs c an be used as SFH 
prob es into the very dist a nt Un i verse ( Lamb & Rc ichartl 
200d; iPorciani fc Madaul l200ll; iBromm fc Loebl |2002|; 



and GRB 090423 , at z = 8.26 jSalvaterra et all 120091 ; 
iTanvir et afll2009ft . 

In principle, the redshift distribution of GRBs can give 
us important clues on the early stages of cosmic history. In 
practice, the connection between GRBs and the underlying 
host galaxy star form ation mode is far from trivial (see e.g. 
iKocevski etafl l|2009l )). thus making the probe value subject 
to uncertainties. 

The purpose of this work is to show how Prin- 
cipal Component Analysis (PCA) can be used to map 
the GRB redshift distribution onto the cosmic SFH 
in a model-independent way. PCA has already been 
applie d to other astrophysical and cosmological con- 



Tot aniet all [20061; iFvnbo et al.l |20( 
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2007; 



Price et al.l 
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potentially to higher redshifts than 
allowed by galaxies alone. For example, GRB 090429B at 
2 = 9.4 (|Cucchiara et al.ll201ll ) is the current record-holder 
object, followed by a z — 8.6 galaxy (|Lehnert et al.ll2010l ) 



texts (iHuterer fc Starkman 2003; Mat uri fc Mignond [20091 : 
iMitra et alj|2010l ; llshida fc de Souzall201ll ). and recognized 
as a useful tool to reconstruct parameters without the in- 
troduction of ad hoc parameterization^]. 



2 THEORETICAL GRB RATE 

We assume that the formation rate of long GRBs 
(duration longer than 2 sec) follows closely the SFH 
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1 Throughout the pa per we adopt a ACD M model with WMAP7 
best fit parameters l ljarosik et al.l l201lh . Q m = 0.267, n A = 
0.734, and H = 71 km s" 1 Mpc _1 . 
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e.g. iTotanil dl997h: ICampisi et all d2010f): IcTardi fc Loebl 



J2000I); ICampisi et all J201ll); IConselice et al.1 (|2005h , 
iBromm fc Loebl (|2006l ): Ide Souza et alJ l|201ll ))~Hence. the 

comoving rate of observable GRBs, *grb is 



i(z) = (n oba /An)f GB .BfbPzp*(z) 



$(L)dlogL, 
(1) 

where £l b s is the field of view of the experiment, /grb is 
the GRB formation efficiency, f b is the beaming factor of the 
burst, P z is the fraction of GRBs with measured redshift and 
Lnm(z) is the luminosity threshold of a given experiment. In 
the following we discuss in detail each of the terms in eq. 
©• 

We define telescope-relate d quantities in eq. ffl after 
Swift specs, i.e. Q b s = 1.4 ijSalvaterra et al.1 120081 ) and 
P z = 0.24 ± 0.060 The overall GRB rate depends on 
ft « 6 2 12, where i s the opening angle of the jet. Accord- 
ing to lGuetta et al.1 (|2005l ) the avera ge value of fb = 0.01 — 
0.02. Using a radio transients survey. [Gal- Yam et al.1 (|2006T l 
placed the upper limit f b < 0.016. We set f b = 0.013 ±0.003 
as a fiducial value. 

The stellar Initial Mass Function (IMF), </f>(m), de- 
termines the fraction of stars massive enough to leave 
a black hole remnant. Current theories indicate that 
the t hreshold mass t o tri gger a GRB is Mgrb = 
25M^ (|Bromm fc Loebl |2006h. Howev er, only a fraction 
Cgrb » 10~ 3 l|Langer fc NormanF 2006) of black holes result- 
ing from supernova explosion actually gives rise to a GRB; 
to be conservative, we considered Cgrb = (1.0±0.5) x 10 -3 . 
Hence, the GRB formation efficiency factor per stellar mass 
is 



/grb = Cgrb 



J^ up m<j>(m)dm 



(2) 



For simplicity, we assume a "standard" Salpeter IMF, 
(m) oc rn,- 235 wit h (M low ,M up ) = (O.lA/ , 1OOM ) 
Schneider et alj|2006t ). 

The number of detectable GRBs depends on the in- 
strument sensitivity and intrinsic isotropic GRB luminosity 
function. F or the latter, we adop t a po wer-law distribution 
function of IWanderman fc Piranl (|2010l ) 



#(L) = 



L < L,, 
L > L». 

1+0.2 



(3) 



with L, = 10 52 b erg s _1 , a L = 0.2+^;^ and /3 L = 1.4±g;|. 
The luminosity threshold is then Ln m = 4-7T c?l Fj^ , where 
d,L is the luminosity distance and Fu m is the bolometric en- 
ergy flux limit of the instrument. In w hat follows we set 
F lim = 1.2 x 10~ 8 erg cm" 2 s" 1 |l1 |2008| ). 

From eq. (JTJ) we can determine the number of observed 
GRBs with redshift in (2, z + dz) over a time interval, At, 
in the observer rest frame: 



dNgm 
dz 



GRB 



At dV 

1 + z dz ' 



(4) 



2 Following IWanderman &: Piranl j2010l ) we consider redshift 
measurements obtained via absorption and photometry only. For 
further details see Sec. l4l 



where dV/dz is the comoving volume element per unit of 
redshift. 

Our observable, i.e. the cumulative number N(z) of 
GRBs up to redshift z, is given by 



N(z) 



dz' 



''-dz' . 



(5) 



So far we have discussed the physical meaning of all terms 
in eq. {1} except the SFH. Since our aim is to build a model 
which is independent from the specific form of p*(z), we 
avoid making hypothesis about this quantity. Instead, our 
intention is to derive p*(z) by using PCA and a (mock or 
real) data set whose data points are {zi, N(zi)} pairs. 



3 PRINCIPAL COMPONENT ANALYSIS 

The main goal of PCA is the dimensionality reduction of 
an initial parameter space through the analysis of its inter- 
nal correlations. Suppose we have a model composed by P 
parameters. If two of them are highly correlated, they are 
actually providing the same information. This means that 
it is possible to rewrite the data in a new parameter space 
consisting of P — 1 terms, wi th minimum l oss of informa- 
tion (for a complete review see IJolliffel (|2002h V This new set 
of parameters are recognized as the principal components 
(PCs), or the eigenvectors of the Fisher information matrix, 
F. 

We postulate that the data set is composed by JV inde- 
pendent observations, each one characterized by a Gaussian 
probability density function, fi\g{xi,fi); Gi, <Tj]. In our no- 
tation, Xi is a measurement of an independent variable; d 
represent the measurements of a quantity G depending on 
Xi\ Gi is the uncertainty associated with the measures, and 
f} is the parameter vector of the theoretical model. In other 
words, we investigate a specific quantity, g, which can be 
written as a function of the parameters Pi. In this context, 
the likelihood function is given by L = nfli /» an d the 
Fisher matrix is defined as 



Fm 



d 2 In L(J3) \ 



(6) 



Brackets in eq. ((6| represent the expectation value. 

We can now diagonalize F, and determine the set of 
its eigenvectors/PCs, e, and eigenvalues, A. Following the 
standard convention, we enumerate e; from the largest to 
the smallest associated eigenvalue. Our ability to determine 
the form of each PC is given by ope; = \ . 

The set e forms a complete base of uncorrelated vectors. 
This allows us to use a subspace of e, eM, to rewrite g as a 
linear combination of all the elements in eM, g rec (x,a). The 
data is then used to find the values of the linear expansion 
coefficients, a (for a detailed discussion see 12011). 

The question of how many PCs should be used in the 
final reconstruction, or how to choose the dimensionality of 
eM, depends on the particular data set analyzed and our 
expectation towards them. To provide an idea of how much 
of the initial information (variance) is included in our plots, 
we shall order them following their cumulative percentage 
of total variance. A reconstruction with the first M PCs 
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encloses a percentage of this value 



t M = 100 



It is important to emphasize that each added PC brings 
its associated uncertainty (<Jpd) into the reconstruction. 
So, although the best-fitted reconstruction converges to the 
"real" function as M increase^, the uncertainty associated 
also raises. As a consequence, the question of how many PCs 
turns into a matter of what percentage of total variance we 
are willing to enclose. 



3.1 Star formation history from GRB distribution 

To specify our method of SFH reconstruction from GRB 
data, let us consider a data set formed by T measure- 
ments of the cumulative number of GRB up to redshift Zi, 
Ndatai = N(zi), and its corresponding uncertainty (<Tj). The 
likelihood is given by 



N datai -N{ Zi ,p) 



Using equation (|6}, the Fisher matrix components are 
1 dN{zi,p) dN{zi,P) 



d0i 



(8) 



(9) 



The Fisher matrix determination is now a matter of calculat- 
ing the derivatives of N(z,/3), for which we use the theoret- 
ical prescriptions of Sec. [2] Aiming at model independence 
and simplicity, we model the SFH as 



(10) 



where j3i are constants, ?ibin is the total number of redshift 
bins, and Ci(z) is a window function which returns 1 if Zi < 
2 ^ Zi+\ and otherwise. Using this description, we may 
write any functional form with resolution limited by our 
computational power. 

The derivatives of N(z) can be computed analytically: 



d0 k 



= H( J{z) + 1 - k) 



J(z) 

} J S k ,iA(zh, zl i+ i)+ 
1=1 



(11) 



+8k,j(z)+iA( z h(z)+i,z)\ , 

where H(x) is a step function which returns if x < and 
1 otherwise, J{z) corresponds to the number of integer bins 
up to redshift z, Si.j is the Kronecker delta function, zU is 
the lower bound of the i-th redshift bin, and 

A / \ ^obsAt r r j-, 
A(Z 1 ,Z 2 ) = JbJGRBfz X 

(12) 



3 This limit holds in ideal situations, where the number of data 
points largely exceeds the number of initial parameters. If the 
data set is small or approximately the same size of the initial 
parameter space, one should also care for overfitting problems. 
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Figure 1. PC A reconstructions of SFH obtained from our mock 
data, using 1 (top-left) to 6 (bottom-right) PCs. The blue-thin 
line corresponds to our fiducial model pfid- The black-thick line 
is the final reconstruction for each case and the red-dashed-thick 
lines corresponds to 2a confidence levels. The inset shows the 
cumulative percentage of total variance, tM- 



From these relations the Fisher matrix can be computed and 
the functional form of the SFH reconstructed through PCA. 



4 RECONSTRUCTION 

Mock data The mock sample is composed of data pairs 
{zi, N(zi)}, distributed in redshift bins of width Az — 0.1, 
where z% represents the middle of each bin. The fiducial 
model used for the SFH is a simple double-e xpon ential func- 
tion, pfld, fitted to numerical results by[Lj (|2008T ). 



logpfld(z) = a + blog (1 + z), 



(13) 



with (a,b) = (-1.70,3.30) for z < 0.993; (-0.727,0.0549) 
for 0.993 < z < 3.80 and (2.35, -4.46) for z > 3.80. We gen- 
erated 500 simulations, each realization with uncertainty in 
the determination of N(z) set to unit and containing 65 
redshift bins. After we generated the mock sample, the in- 
formation about /Ofid(z) is discarded. Our goal from now on 
is to re-obtain the functional form of eq. (|13jl from PCA. 

The Fisher matrix is obtained as described above, as- 
suming an observing time A£ms = 1 yr. As the mock sample 
purpose is to test the procedure under an ideal scenario, we 
did not include uncertainties in the parameters of eq. (flj. 
Having obtained the PCs, we can rewrite the SFH as 



prec\Z J 



M 

1=1 



(14) 



where ctj and p c are constants to be determined and M is the 
number of PCs we choose to use in the reconstruction. The 
simulated data points are then used to find the appropriate 
values for the parameters aj and p c as those that minimize 
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the expression 

X 2 (a)ocX: ( ^ data -y^ e ' Q))2 , (15) 

i — 1 8 

where <7j = 1 for all redshift bins. The reconstructions ob- 
tained using 1 to 6 PCs are shown in Fig. [T] The uncertainty 
in the final reconstruction was calculated by a quadrature 
sum that includes the parameters apd and the uncertainty 
in the determination of parameters a(a ai ) and p c (o"p c ). 

From Fig. Q] we can appreciate the success of the proce- 
dure in reconstructing the underlying unknown SFH in an 
ideal scenario, with increasing agreement as the number of 
PCs raises. Confidence levels also become wider as M in- 
creases, with the only exception of the reconstruction with 
1 PC. Since a Pa dominates the errors due to the limited free- 
dom to fit the second peak of the fiducial model with only 1 
PC. With 2 PCs fitting the second peak becomes easier, and 
as a consequence, the magnitude of a Pc decreases to levels 
below those of ope; • 

Swift data After validating PCA reconstruction under 
ideal conditions, we turn to the use of currently available 
Swift data, and compare these results with independent 
measurements of SFH from the literature. First, we need 
to properly choose our data set. Since only GRBs with mea- 
sured redshifts can be used in our analysis, the question of 
how the redshift measurements were obtained must be ex- 
amined carefully. 

GRBs redshifts are generally obtained from optical af- 
terglow spectra using absorption lines or photometry, or 
from the spect rum of the host galaxy usin g emission lines. 
As pointed bv lWanderman fc Piranl (|2010l ). different meth- 
ods yield different redshift distributions: a visual inspection 
of Fig. [2] illustrates this point. Most noticeably, the GRB 
redshift distribution determined from their hosts lacks very 
high-z events. Moreover, emission (and to a lesser extent, ab- 
sorption) lines are susceptible to a selection effe ct known as 
the "redshift desert " in the range 1.1 < z < 2.1 |Fiore et al.l 
l2007l ; ICowar"dll2009l). Additional bia s sources are preliminary 
discussed bv lMalesani et al.l (|2009f ). To avoid systematic er- 
rors affecting the overall redshift distribution, our data sam- 
ple is composed by 120 Swift GRBs with redshift determined 
from absorption lines and photometry (gray region in Fig. 
[21 top panel). 

The next step is to choose the appropriate redshift bin 
width. In principle, the quality of the reconstruction should 
increase with the number of bins. However, as GRB are dis- 
crete events, if we pick a bin width based on the available 
data (for example, in such a way that each bin has at least 
one GRB), the bins will be too wide (~ 1). In this case, 
the assumption that the SFH is constant inside the bin will 
not hold, leading to reconstructions with bad resolution. 
To overcome this limitation we performed a gaussian ker- 
nel fi10 to the data (black line in Fig. [2l top panel). Now we 
have a continuous probability distribution function (PDF) 
for dN/dz, which follows the real data distribution and al- 

4 A non-paramctric estimate of the PDF obtained from a linearly 
interpolated version of ^ 2~Z?=i ^ ( X ~h* J f° r a kernel k(x), bin 
width h and a total of n bins. 
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Figure 2. Top: Histogram showing the measured redshift distri- 
bution of 120 GRBs detected by Swift from 2005 to 2010, divided 
by redshift measurement methods: absorption (red-dashed), pho- 
tometry (green-dotted) and emission (bluc-dot-dashed). The gray 
region corresponds to the redshift distribution of data we used 
(absorption + photometry). The solid black line shows the fit 
to the used data distribution. Bottom: The cumulative distribu- 
tions constructed from Swift data (gray stars) and from our fitted 
distribution function (black points). 

lows us to set the bin width as small as required. We kept 
Az = 0.1 and use the PDF to calculate the cumulative num- 
ber of observed GRBs in each bin. The comparison between 
the real data cumulative distribution and the one calculated 
via the fitted PDF are shown in Fig. [2j bottom panel. 

The Fisher matrix is calculated using Atsw = 6 yr of 
observation time. The parameters Oi were obtained by sum- 
ming in quadrature the uncertainties in the quantities in- 
volved in eq. fl}. Fig. [3] shows the first 2 PCs and the corre- 
sponding reconstruction using both of them, which already 
encloses more than 97% of total variance. In the lower panel, 
the points correspond to completely independent measure- 
ments from the literature. These data points are shown only 
for comparison purposes and have not been used in our cal- 
culations. 



5 DISCUSSION 

We have proposed the use of PCA as a powerful tool to re- 
construct the cosmic star formation history exploiting the 
measured gamma-ray burst redshift distribution. The pro- 
cedure was successfully validated using synthetic data and 
next applied to actual Swift data (Fig. 3). 

It is important to highlight that the approach is com- 
pletely independent of the initial choice of the theoretical 
model parameter vector, /?. This has the obvious advantage 
of avoiding any a priori hypothesis on the SFH, p*(z). How- 
ever, the degeneracy between p»(z) and any other factor 
we are failing to take into account cannot be removed, i.e. 
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Figure 3. Top: First (red solid) and second (blue dashed) PCs 
from Swift. Bottom: PCA reconstruction from Swift data using 
2 PCs, compared with independent SFH determinations (light 
points, not used in our calculations). The black solid line is the 
PCA best-fit reconstruction using 2 PCs; the red dashed lines cor- 
respond to 2a confidence levels. The inset shows the cumulative 
percentage of total variance, tj^f • 



the reconstructed SFH contains also the behavior of all the 
agents influencing the determination of N(z) and not in- 
cluded in the model. 

For example, lLanger fe Normanl (|200fj ); 

IWooslev fe Hegerl (|2006l ) have argued that GRB progenitors 
will have a low metallicity. Such an effect would be a con- 
sequence of the mass and angular momentum loss caused 
by winds in high-metallicity stars. This would prevent such 
stars of becoming GRBs and consequently, change their ex- 
pected redshift distrib utio n dSalvaterra fe Chincarinil 120071 ; 
ISalvaterra et all 120071 ; O 120081 ). We implicitly considered 
that this and other such effects will span within the error 
bars in our analysis. 

In spite of the remaining uncertainties, which are prob- 
ably less severe than those affecting other methods aimed at 
recovering the high-z tail of the SFH, there are robust indica- 
tions that we can gather from the analysis of our results. The 
first is that the combination of GRB data and PCA suggest 
that the level of star formation activity at z sa 9.4 could have 
been already as high as the present-day one (« O.OIM0 yr _1 
Mpc~ :! ). This is a factor 3-5 times higher than deduced from 
high-z galaxy searches thr ough drop-out techniqu es, simi- 
larly to the trend found by lYonetoku et aU (|2004l ). If true, 
this might alleviate the long-standing problem of a photon- 
starving reionization; it might also indicate that galaxies ac- 
counting for most of the star formation activity at high red- 
shift go undetected by even the most deep searches. Finally 
it is worth noticing that a sustained high-z star formation 
activity is consistent with predictions of reionization mod- 
els that simultaneously match a number of observable ex- 
perimental constraints as the Gunn-Peterson effect, Thom- 



son free-e lectron optical depth, Lyman L imit Systems statis- 
tics et c. (|Choudhurv fc Ferraral l|200fj ). iBolton fc HaehneltJ 
(2007)). Given the expected large growth of GRB detec- 
tions from the next generation of instruments, the method 
proposed here promises to become one of the most suitable 
and reliable tools to constrain the star formation activity in 
the young Universe. 
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